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DIRECT  BLIND  DECONVOLUTION  II. 
SUBSTITUTE  IMAGES  AND  THE  BEAK  METHOD 


ALFRED  S.  CARASSO* 

Abstract.  The  BEAK  method  is  an  FFT-based  direct,  blind  deconvolution  technique  previously 
introduced  by  the  author,  and  applied  to  a limited  but  significant  class  of  blurs  that  can  be  expressed 
as  convolutions  of  2-D  symmetric  Levy  probability  density  functions.  This  class  includes  and  general- 
izes Gaussian  and  Lorentzian  distributions,  but  does  not  include  defocus  blurs.  The  method  requires 
a-priori  information  on  the  Fourier  transform  /,  (£,  '/)  of  the  unknown  exact  image  fe(x,y),  namely, 
the  gross  behavior  of  log  |/e  (£,  h)l  along  a single  line  through  the  origin  in  the  (£,?/)  plane.  The 
present  paper  significantly  extends  the  applicability  of  the  BEAK  method.  It,  is  shown  that  images 
of  similar  objects  often  display  approximately  equal  gross  behavior,  and  that  gross  behavior  in  such 
substitute  images  can  be  used  successfully  in  numerous  practical  contexts.  Next,  using  substitute 
images,  a variant  of  the  BEAK  method  is  developed  that  can  handle  defocus  blurs.  The  paper  is 
illustrated  with  several  examples  of  blind  deconvolution  of  512  x 512  images  in  the  presence  of  noise, 
and  includes  a detailed  discussion  of  an  example  where  the  BEAK  method  fails. 

Key  words,  image  deblurring,  blind  deconvolution,  direct  methods,  Levy  density  functions, 
defocusing,  substitute  images,  BEAK  method,  SECB  method. 

AMS  subject  classifications.  35R25,  35B60,  60E07,  68U10. 


1.  Introduction.  This  paper  is  a sequel  to  [6].  As  was  the  case  there,  the 
procedure  described  below  will  generally  not  be  useful  for  severely  blurred  images 
at  high  levels  of  noise,  nor  for  images  degraded  by  arbitrary  or  unknown  processes. 
Rather,  the  method  is  limited  to  a narrow  class  of  deblurring  problems  involving 
restricted  types  of  blur  at  low  to  moderate  intensities,  in  the  presence  of  low  levels  of 
noise.  With  the  type  of  blur  assumed  known , these  conditions  enable  identification  of 
the  parameters  in  the  system  point  spread  function.  The  discussion  below  includes 
an  example  where  the  method  fails.  Other  examples  of  failure  are  easily  found.  All 
images  in  this  paper  are  of  size  512  x 512.  Unless  otherwise  indicated,  these  images 
are  quantized  at  8-bit,s  per  pixel,  i.e.,  each  pixel  value  is  an  integer  lying  between  0 
and  255.  Such  quantization  introduces  8-bit  rounding  noise  which  plays  a significant 
role.  Additional  noise  processes  are  applied  in  some  cases. 

Blind  deconvolution  seeks  to  deblur  an  image  without  knowing  the  point  spread 
function  describing  the  blur.  In  [6],  two  methods  were  developed  for  direct  (i.e.,  non- 
iterative) blind  deconvolution,  the  BEAK  method  and  the  APEX  method.  These 
are  Fourier  domain  techniques  for  detecting  the  signature  of  the  system  point  spread 
function  from  1-D  analysis  of  the  blurred  image.  A separate  direct  Fourier  domain 
deblurring  technique,  the  SECB  method  [4],  uses  this  detected  point  spread  function  to 
deblur  the  image.  Direct  blind  deconvolution  of  512  x 512  images  can  be  accomplished 
in  minutes  on  current  desktop  workstations.  Other  approaches  to  that  problem  are 
iterative  in  nature,  but  the  iterative  process  is  not  always  well-behaved  [6].  When  that 
process  is  stable,  several  thousand  iterations  may  be  necessary  to  achieve  convergence 
[11],  [15].  The  discussion  in  [4]  and  [6]  is  focused  exclusively  on  so-called  class  G point 
spread  functions.  This  is  a significant  but  limited  class  of  blurs  that  can  be  expressed 
as  convolutions  of  2-D  symmetric  Levy  ‘stable’  probability  density  functions.  This 
class  includes  and  generalizes  Gaussian  and  Lorentzian  distributions,  but  does  not 
include  defocus  blurs.  The  BEAK  detection  method  uses  a novel  type  of  a-priori 
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information  about,  the  exact  sharp  image  fe(x,y),  namely,  the  gross  behavior  of  its 
Fourier  transform  /e(£,  //)  along  a single  line  through  the  origin  in  the  (£,77)  plane. 
Such  gross  behavior  is  represented  by  a smooth  function  that  provides  a least  squares 
fit  to  the  highly  oscillatory  trace  of  log  |/e(£,7/)|  along  that  line.  By  pre-rotating 
the  image  if  necessary,  one  may  assume  that  when  such  gross  behavior  is  known,  it 
is  known  along  the  horizontal  axis  (£,0).  A large  class  of  sharp  images,  the  class  W, 
is  exhibited  in  [6]  having  the  property  that  gross  behavior  can  be  summarized  by  the 
two  positive  numbers  a,  b in  the  expression 

(1)  log  |/e*(£,  0)|  ~ ~a  ICI2 * * * 6,  «,  b > 0,  fe{x,y)  G W. 


where 

(2)  |//(d0)|  = \fe  (£,  0)  | / f e (0,  0) . 

The  second  detection  method  discussed  in  [6],  the  APEX  method,  does  not  require 
prior  information  on  fe(x,y)  but  assumes  that  image  to  be  a recognizable  object. 
Several  (fast,)  interactive  trials  are  necessary  before  locating  a suitable  point  spread 
function  using  that  method. 

In  the  present  paper,  the  APEX  method  is  not  discussed  further.  Rather,  the 
range  of  applicability  of  the  BEAK  method  is  extended  in  two  major  ways.  First, 
we  observe  that  gross  behavior  in  class  W images  provides  only  superficial  prior 
information  about  /,,  (:r. ;//),  and  that  images  of  similar  objects  are  often  found  to 
display  approximately  equal  gross  behavior.  Therefore,  in  using  the  BEAK  method 
for  identifying  system  point  spread  functions,  it  is  generally  not  necessary  to  know  the 
gross  behavior  in  the  original  image  fe(x,  y),  but  only  that  in  /s(.x,  y),  an  appropriate 
sharp  image  of  a.  similar  object.  There  are  numerous  practical  contexts  where  such 
similar  images  are  available.  By  substituting  the  gross  behavior  in  fs(x,y)  for  that 
in  the  unknown  fe(x,y ),  the  BEAK  method  becomes  applicable  in  a wide  variety  of 
situations.  This  is  illustrated  with  examples  of  class  W images  blurred  by  class  G 
point,  spread  functions. 

A second  aim  of  this  paper  is  to  develop  a detection  method  for  defocus  blurs. 
As  will  be  seen  below,  deconvolution  of  defocus  blurs  is  strikingly  different  from  class 
G blurs.  Using  substitute  images  as  described  above,  a variant  of  the  BEAK  method 
is  presented  that  can  approximately  identify  defocus  point  spread  functions,  provided 
the  defocus  is  not  too  severe.  This  detected  psf  can  then  be  used  to  deblur  the  image. 
Since  the  SECB  deblurring  procedure  used  in  [6]  is  restricted  to  class  G blurs,  a, 
modification  of  that  procedure  is  necessary  for  defocus  and  other  shift -invariant  blurs 
not,  in  G.  This  modification  is  discussed  in  the  Appendix. 

2.  Gross  behavior  in  class  W images.  A-priori  information  is  paramount,  in 

the  solution  of  ill-posed  inverse  problems.  Such  information  becomes  ever  more  critical 

in  the  case  of  blind  deconvolution,  where  severe  non-uniqueness  is  compounded  with 
discontinuous  data  dependence.  I11  iterative  blind  deconvolution  algorithms,  positivity 
and  support,  constraints  011  the  convolution  components  are  generally  imposed  in  order 

to  reduce  the  multiplicity  of  solutions.  However,  such  constraints  are  not,  always 
effective  in  avoiding  traps  at  local  minima  (stagnation  points),  or  divergence  of  the 

iterative  procedure.  There  is  therefore  considerable  interest  in  finding  types  of  prior 
information  that  are  widely  applicable,  and  lead  to  reliable  algorithms.  Gross  behavior 
in  Fourier  space  is  one  example  of  useful  prior  knowledge. 
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Fig.  1.  log  |/e  (£,0)|  ^ —a  |£r  and  rough  equivalence  of  gross  behavior  in  two  different. 


spacecraft  images.  A.  Mariner  5 has  a = 2.81,  b 
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SHARP  US'S  EISENHOWER  IMAGE 


DISCRETE  FREQUENCY  {INTEGER) 
LEAST  SQUARES  PIT  ( SOLID  LINE | 


0.190  . B.  Mariner  10  has  a = 2.98,  b = 0.183. 
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SHARP  USS  KITTY  HAWK  IMAGE 
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DISCRETE  FREQUENCY  ( INTEGER ) 

LEAST  SQUARES  PIT  (SOLID  LINE I 


Fig.  2.  log  |/e  (/,0)|  ~ —a  |/|b  an(l  Tough,  equivalence  of  gross  behavior  in  two  different 
aircraft  carrier  images.  A.  USS  Eisenhower  has  a = 2.91,  b — 0.157  . B.  USS  Kittyhawk  has 
a = 3.04,  b - 0.158. 


In  the  simplest  case,  image  deblurring  is  associated  with  the  solution  of  two- 
dimensional  convolution  equations 

(3)  Hf  = / h(x  - u,  y - v)f{u,  v)dudv  = h{x.  y)  ® f(x.  y)  = g{x,y ), 

Jr- 

where  g(x.y)  is  the  recorded  blurred  image,  f(x,y)  is  the  desired  unblurred  image, 
h(x,y)  is  the  blurring  kernel  or  point  spread  function  (psf)  of  the  imaging  process, 
and  (g>  denotes  convolution.  The  psf  h(x,y)  represents  the  cumulative  effects  of  all 
distortions  caused  by  the  media  through  which  signals  propagate,  as  well  as  all  optical 
and  electronic  aberrations  produced  by  imperfect  sensing  and  recording  equipment. 
It  is  assumed  that  h(x,y)  is  such  that  the  linear  problem  Hf  = g has  at  most  one 
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SHARP  UARIIYN  MONROE  IMAGE 


SHARP  INGRID  BERGMAN  IMAGE 


FiG-  3.  log  |/e  (£,0)|  fts  —a  |£|6  and  rough  equivalence  of  gross  behavior  in  two  different  face 
images.  A.  Marilyn  Monroe  has  a — 3.57,  b = 0.158  . B.  Ingrid  Bergman  has  a = 3.11,  b = 0.183. 
Solid  curve  in  A deviates  only  slightly  from  solid  curve  in  B,  despite  differences  in  corresponding 
values  for  a and  b. 


solution.  This  is  the  case  for  class  G point  spread  functions  defined  in  Fourier  space 

by 

(4)  h(£,Ti)  = e-'E>'=iai(t2+r,3)0',  a,>  0,  0 < fa  <1. 

The  function  h(£,r])  is  known  as  the  optical  transfer  function  (otf)  of  the  imaging 
process.  The  blurred  image  g(x,y ) includes  noise,  which  is  viewed  as  a separate 
additional  degradation, 

(5)  g(x,y)  = ge{x,y)  + n(x,y), 

where  ge(x,y)  is  the  blurred  image  that  would  have  been  recorded  in  the  absence  of 
noise,  and  n{x,y ) represents  the  cumulative  effects  of  all  noise  processes  and  other 
errors  affecting  final  acquisition  of  the  digitized  array  g(x,  y).  This  includes  nonlinear 
noise  processes  where  n{x,y)  may  be  a function  of  f(x,y).  Both  ge(x1y)  and  n(x,y) 
are  unknown,  but  n(x,y)  may  be  presumed  small.  The  unique  solution  of  (3)  when 
the  right  hand  side  is  ge{x,y ),  is  the  exact  sharp  image  denoted  by  fe(x,y).  Thus 

(6)  h(x,  y)  ® fe(x,y)  = ge(x,  y). 

Since  fe(x,y)  > 0 

(7)  \fe{^v)\<  I fe{x,y)dxdy  = /e(0,0)  = a > 0. 

JR 2 

Also,  since  h(x.y)  is  a probability  density, 

(8)  ge(0,  0)  = / ge(x,  y)dxdy  = [ fe(x,y)dxdy  = fe{0, 0)  = a > 0. 

jr2  Jr 2 

Using  <7  as  a normalizing  constant,  we  may  normalize  Fourier  transform  quantities 
q(£,i  V)  by  dividing  by  a.  Let 


(9) 
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Fig.  4.  log|/e  (£,0)|  ~ —a  |£|b  and  gross  behavior  in  three  different  MRI  brain  images. 
A.  Brain  1 has  a = 3.12,  b = 0.155.  B.  Brain.2  has  a = 3.28,  b = 0.154.  C.  BrainS  has 
a — 4.03,  b = 0.135.  D.  Gross  behavior  roughly  equivalent  in  Brainl  and  Bram2 , but  markedly 
different  in  Brainl  and  BrainS.  Use  of  BrainS  as  substitute  for  Brainl  in  BEAK  method  leads  to 
significantly  better  results  than  does  use  of  BrainS.  See  Fig.  8. 


denote  the  normalized  quantity.  The  function  | fe  (£,r/)|  is  highly  oscillatory,  and 
0 < |/e  | < 1.  Since  fe(x,y)  is  real,  its  Fourier  transform  is  conjugate  symmetric. 
Therefore,  the  function  |/e  (£,  ?7)|  is  symmetric  about  the  origin  on  any  line  through 
the  origin,  77  = £tan$,  in  the  (£,17)  plane. 

As  previously  noted,  all  images  in  this  paper  are  of  size  512  x 512.  For  each  sharp 
image  /e(x,  y),  the  discrete  Fourier  transform  is  a 512  x 512  array  of  complex  numbers, 
which  we  again  denote  by  /e(£,  77)  for  simplicity.  The  ‘frequencies’  £,  ?/  are  now  integers 
lying  between  —256  and  256,  and  the  zero  frequency  is  at  the  center  of  the  transform 
array.  This  ordering  is  achieved  by  pre-multiplving  fe(x.y)  by  ( — l)J+y.  The  discrete 
transforms  fe  (£,0),  and  fe  (0,7 7)  are  immediately  available.  Image  rotation  may  be 
used  to  obtain  discrete  transforms  along  other  directions.  All  1-D  Fourier  domain 
plots  shown  in  this  paper  are  taken  along  the  axis  77  = 0 in  the  (£,77)  plane.  In  these 
plots,  the  zero  frequency  is  at  the  center  of  the  horizontal  axis  and  the  graphs  are 
symmetric  about  the  vertical  line  £ = 0. 
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The  class  W of  well-behaved  sharp  images  f,(x,y)  introduced  in  [C]  may  be  loosely 
characterized  as  follows: 


* 

1.  | f,  (£,//)  | has  at  most  isolated  zeroes  in  the  (£,77)  plane. 


2.  Neglecting  isolated  singularities,  the  global  behavior  of  log  | fe  (£,?/) | on  any  line 
i/  — £ tan 0 is  roughly  monotone  decreasing  with  increasing  r = (£2  + 772)1/2;  i.e.,  on 


the  ray  re’0,  a least  squares  tit.  to  log \fe  (£,t/)|  with  an  appropriate  monotone  de- 
creasing  function,  provides  a fair  representation  of  gross  data  behavior  as  r increases. 

3.  While  the  rate  of  decay  may  vary  between  rays,  this  decay  is  relatively  slow,  i.e.,  of 
a general  order  of  magnitude  comparable  to  that  found  in  the  sharp  Mariner  images 
in  Fig  1. 

4.  The  gross  behavior  of  log  \ fe  (£,  7/)|  along  the  ray  re’0  is  defined  to  be  the  function 

v(r)  = —a(6)  , a,h  > 0,  that  best  fits  log|/P  (£,£tan#)|  in  the  least  squares 


sense. 


Nine  examples  of  class  W images  are  displayed  in  Figs.  1 through  4,  together  with 
their  gross  behavior  along  the  axis  77  = 0 in  the  (£,77)  plane.  Several  other  examples 
are  given  in  [6].  What  is  of  interest  here,  is  the  fact  that  ‘similar’  objects  appear  to 
have  roughly  equivalent  gross  behavior.  In  Figs.  1 and  2,  the  values  for  a and  b are 
approximately  equal  in  the  two  spacecraft  and  in  the  two  carrier  images.  Hence,  the 
corresponding  gross  behavior  traces  are  roughly  equivalent.  In  Fig.  3,  the  values  of  a 
and  b for  Marilyn  Monroe  and  Ingrid  Bergman  are  substantially  different.  However, 
the  corresponding  gross  behavior  traces  almost  coincide  in  that  case.  In  Fig.  4,  gross 
behavior  in  Brain3  is  markedly  different  from  that  in  Brain  1 and  Brain2.  This  shows 
that  similar  objects  need  not  always  have  equivalent  gross  behavior.  In  practice,  it 
may  be  necessary  to  sample  several  similar  objects  before  arriving  at  a good  substitute 
image  for  use  in  the  BEAK  method. 

3.  The  BEAK  method  for  class  G psfs.  As  discussed  in  [6],  this  is  a Fourier 
domain  technique  for  detecting  class  G point  spread  functions  acting  on  class  W 
images.  The  method  uses  1-D  Fourier  analysis  of  the  blurred  image  data  g(x,y ) 
in  (3),  and  requires  prior  knowledge  of  gross  behavior  in  the  unknown  sharp  image 
fe{x,y ) in  (6).  The  detected  optical  transfer  function  is  then  used  to  solve  the  ill-posed 
problem  (3).  This  is  accomplished  using  another  direct  Fourier  domain  procedure,  the 
SECB  method  [2,  3,  4], 

The  BEAK  method  is  based  on  the  following  observations.  In  the  basic  relation 


(10) 


y(x,  y)  = h{x.  y)  0 fe(x,y ) + n(x,y ), 


we  may  safely  assume  that  the  noise  n(x,y ) satisfies 


(ID 


so  that, 


(12) 


n*(£,t  7)1  « 1- 
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DISCRETE  FREQUENCY  ( INTEC-ER I DISCRETE  FREQUENCY  ( INTECER ) 

LEAST  SQUARES  FIT  iSOUD  LINE I 


Fig.  5.  Application  of  BEAK  method  to  8-bit  blurred  Eisenhower  image,  using  gross  behavior 
from  sharp  Kittyhawk  image.  A.  8-bit  blurred  Eisenhower  image  g(x,y)  obtained  by  convolution 
of  sharp  Eisenhower  image  fe(x,y)  with  'long  exposure  turbulence ’ Levy  psf  with  a = 0.003,  0 = 
5/6.  Trace  of  log  |<;*  (f , 0)|  above  8-bit  noise  level  only  on  small  interval  —60  < £ < 60.  B. 
log  |/s  (e,  0)|  ss  —3.04  |£|0158  in  sharp  Kittyhawk  image  will  be  used  to  obtain  BEAK  fit  to  blurred 
Eisenhower  trace. 


Consider  the  case  where  the  otf  is  a pure  Levy  density  /)(£,  ?/)  = e q(C+t)'  Since 

g = ge+n 

(13)  log|r)*(C'd)l  = log  | c“Q(€‘+'r)'3/e*(gO/)  +n*{£,ri) |. 

Let  ft  = {(C??)  | £2  + if  < cu2}  be  a neighborhood  of  the  origin  where 

(14)  e-a^+^yi\ft*(^,))\>>\n*(^V)\. 

Such  an  ft  exists  since  (14)  is  true  for  £ = r\  — 0 in  view  of  (12).  The  radius  lo  > 0 of 
ft  decreases  as  a and  n increase.  For  E 11  we  have 

(15)  log  l.(7*(C  V)\  ~ -0'(^2  +V2)13  + log|/e*(4,//)|. 

Hence,  for  |£|  < lo, 

(16)  log  I//*  (C  0)|  « -a|^  + log|//(C0)|. 

The  idea  in  [G]  is  to  replace  the  unknown  log|/e  (£,  0)|  in  (1G)  by  its  gross 

behavior,  i.e.,  its  least  squares  approximation  v(£)  — — a |£|fc,  which  is  assumed 

known.  Extensive  numerical  experiments  with  a wide  assortment  of  images  and  class 

G psfs,  indicate  this  to  be  a successful  strategy.  In  the  present  paper,  we  replace  the 
~ * , 
unknown  logj/e  (£, 0)|  in  (1G)  by  the  least  squares  approximation  vs(0  — — o |£|" 

to  log  |/s  (£,  0)|,  where  fs  {x,  y)  is  a substitute  image , i.e.,  an  image  of  a similar  object 

that  is  expected  to  have  roughly  equivalent  gross  behavior.  Several  candidate  images 

may  need  to  be  tried  before  achieving  optimal  results.  The  detection  procedure  is 

then  the  following.  With  a and  b given  a-priori,  find  positive  numbers  a,  $, 

so  that  the  function  u(£)  = — a |£|2,2  — a |£|b  best  fits  log  ).(/*(£.  0)|  on  |£|  < lo. 
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Fig.  6.  Blind  deblurring  of  aircraft  carrier  image  using  BEAK  psf  detection  method.  A.  BEAK 
fit  to  8-bit  blurred  Eisenhower  trace  on  —60  < £ < 60,  using  gross  behavior  in  sharp  Kittyhawk  trace. 
Fit  returns  detected  parameters  a = 0.000931,  (3  = 0.97,  differing  from  true  values  a = 0.003,  (5  = 
5/6.  B.  Detected  psf  (dashed)  roughly  approximates  true  psf  (solid).  C.  Blurred  Eisenhower  image. 
D.  SECB  deblurring  of  image  C using  s = 0.001,  K — 1.27,  and  detected  psf  parameters  a = 
0.000931,  (3  = 0.97.  Good  restoration  achieved  even  though  detected  psf  not  equivalent,  to  true  psf. 
Ringing  artifacts  visible  against  light  background  near  image  top. 


This  may  be  accomplished  interactively  using  nonlinear  least  squares  algorithms  in 
DATAPLOT  [10].  The  returned  values  for  a and  /3  are  subsequently  used  for  h(£,  r/) 
in  the  SECB  deblurring  procedure. 

For  more  general  class  G otfs  where  h(£,r))  = e_i:^iQhS'  + 'i"),3!  ^ we  again  seek  the 
best  fit  to  log|$*(f,0)|  on  |£|  < u,  with  a function  u(f)  — — a |£|2^  — a |£|6.  Here, 
the  returned  values  for  a and  (3  may  be  considered  average  values  for  the  a, , and 
are  expected  to  generate  a pure  Levy  density  that  well-approximates  the  composite 
psf. 

We  illustrate  this  procedure  in  Figs.  5 and  6.  In  Fig.  5A,  the  sharp  USS  Eisen- 
hower image  shown  in  Fig.  2A,  was  artificially  blurred  by  convolution  with  a pure 
Levy  density  with  a = 0.003  and  /3  = 5/0.  This  simulates  long  exposure  imaging 
in  the  presence  of  turbulence  [9].  The  result  of  that  numerical  convolution  was  then 
truncated  to  8-bits.  The  trace  of  log  0)|  in  Fig.  5A  lies  above  8-bit  noise  level 
only  on  the  small  interval  —60  < £ < 60,  and  is  obviously  quite  different  from  that 
in  the  original  image  in  Fig.  ‘2A.  The  sharp  LISS  Kittyhawk  image  in  Fig.  5B  will 
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Fig.  7.  Application  of  BEAK  method  to  8-bit.  Brainl  im.age  blurred  with  Levy  density  where 
a = 0.05  and  (3  = 0.6.  Blurred  image  trace  lies  above  8-bit  noise  level  on  interval  —50  < £ < 50.  A. 
BEAK  fit  using  gross  behavior  in  sharp  Brain2  image  returns  detected  parameters  a — 0.0325,  0 — 
0.643,  differing  from  true  values.  B.  BEAK  fit  using  gross  behavior  in  sharp  BrainS  image  develops 
more  curvature,  and  returns  a = 0.00316,  0 = 0.921,  markedly  different  from  true  values.  C. 
Detected  psf  in  A (dotted  curve)  is  much  closer  to  true  psf  (solid  curve),  than  is  detected  psf  in  B 
(dashed  curve). 


Fig.  8.  Blind  deblurring  of  Brainl  image  using  two  distinct  BEAK  detected  psfs.  A.  8-bit 
blurred  Brainl  image.  B.  SECB  deblurring  of  image  A using  s = 0.001,  A = 1.27,  and  Brain2 
detected  psf  parameters  a = 0.0325,  0 = 0.643.  Good  restoration  achieved  even  though  detected  psf 
not  equivalent  to  true  psf.  C.  SECB  deblurring  of  image  A using  s = 0.001.  A = 1.27,  and  Brain3 
detected  psf  parameters  a = 0.00316,  0 = 0.921,  results  in  comparatively  weak  restoration. 


now  be  used  as  the  substitute  image  in  the  BEAK  method.  Using  a = 3.04  and 
b = 0.158.  we  best  fit  log  |g*(£,  0)|  in  Fig.  5A  on  the  interval  |£|  < 00,  with  the 
function  tt(£)  = —a  — a |£|\  using  the  fit  command  in  DATAPLOT  [10].  That 
fit  is  indicated  by  the  solid  curve  in  Fig.  6A,  and  the  beak  or  mb  in  that  curve  near 
£ = 0 is  a characteristic  feature  of  this  detection  procedure.  The  fit  returns  detected 
parameters  d = 0.000931,  (5  = 0.97,  that  differ  substantially  from  the  correct  values 
a — 0.003,  (3  = 5/6.  However,  as  indicated  in  Fig.  6B,  the  detected  psf  (dashed 
curve)  roughly  approximates  the  true  psf  (solid  curve).  Using  this  detected  psf  in  the 
SECB  method  with  s — 0.001  and  I\  = 1.27,  produces  a quite  useful  restoration,  as 
shown  in  Fig.  6D.  While  ringing  artifacts  are  visible  in  the  deblurred  image  against 
the  light  sky  background  near  the  top  of  the  image,  considerable  sharpening  of  the 
planes  on  deck,  as  well  as  the  'island’  structure  containing  the  bridge  and  mast,  has 
been  achieved. 

The  next  example  illustrates  the  importance  of  locating  a good  substitute  image. 
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Tlu'  sharp  MRI  Brain  1 image  in  Fig.  4A  was  blurred  synthetically  by  convolution 
with  a pure  Levy  density  with  a — 0.05  and  [3  = 0.0,  and  the  numerical  convolution 
was  rounded  to  8-bits.  This  time,  the  trace  of  log  |<7*(£,  0)|  lies  above  8-bit  noise  level 
on  tlu-  interval  —50  < £ < 50.  In  Fig.  7A,  the  BEAK  method  is  applied  on  |£|  < 50 
using  Brain‘2  in  Fig.  4B  as  a substitute  image.  This  produces  detected  parameters 
o = 0.0325,  (3  = 0.643.  In  Fig.  7B,  the  BEAK  method  is  applied  using  Brain3  in 
Fig.  1C  as  the  substitute  image.  The  resulting  tit  develops  more  curvature  than  in 
Fig.  7A,  and  returns  detected  parameters  a = 0.00310,  / 3 = 0.921.  These  values  are 
markedly  different  from  the  true  values  a = 0.05,  (3  — 0.0.  As  shown  in  Fig.  7C,  and 
as  may  be  expected  from  Fig.  4D.  the  detected  psf  using  Brain2  is  much  closer  to 
the  true  psf,  than  is  the  detected  psf  using  Brain3.  In  Fig.  8,  we  see  that  the  Brain2 
detected  psf  produces  a quite  good  restoration,  whereas  the  narrower  Brain3  detected 
psf  results  in  a comparatively  weak  reconstruction. 

Numerous  other  experiments  with  class  G psfs  using  the  BEAK  method  with 
substitute  images,  confirm  the  pattern  of  behavior  typified  by  the  above  two  examples. 
It  should  be  noted  that  in  addition  to  locating  a good  substitute  image,  successful 
detection  requires  having  a sufficiently  wide  interval  about  £ = 0 wherein  the  trace 
of  log|r/*(C0)|  lies  above  noise  level.  This  was  the  case  in  Figs.  0A  and  7A.  Severe 
blurring  caused  by  large  values  of  a when  j3  is  near  unity,  and/or  high  levels  of  noise, 
considerably  reduce  the  width  of  that  interval  and  can  preclude  useful  detection. 

4.  Uniform  defocus  blur.  With  class  G blurs,  one  can  often  deblur  the  image 
quite  satisfactorily  with  a psf  that  is  only  a rough  approximation  to  the  true  psf.  This 
was  the  case  in  the  USS  Eisenhower  image  in  Fig.  6D,  the  MRI  brain  image  in  Fig. 
8B,  as  well  as  in  several  other  examples  in  [6].  We  shall  hud  defocus  blurs  to  be  less 
forgiving. 

Defocus  blurs  are  discussed  in  [14].  Important  early  work  on  blind  deconvolution 
of  defocus  blurs  can  be  found  in  [1],  Additional  references  may  be  found  in  [12].  If 
R > 0 is  the  radius  of  the  ‘circle  of  confusion’,  the  psf  for  uniform  defocus  blur  is 
given  by 


f ( 7tT?2  ) 1 , x2  + y'2  < R2, 

(17)  h(x,y)=< 

{ 0,  x2  + y2  > R2 . 

This  has  a Fourier  transform  given  by  the  ‘sombrero  function’  [8,  p.  72] 

(18)  hR,V)  = 2 MR0)/(Rd),  0 = y/£2  + v 2, 

where  J \{x)  is  the  Bessel  function  of  the  first  kind  of  order  1.  A 1-D  cross  section 
of  the  sombrero  function  when  R = 0.08  is  shown  in  Fig.  9A  and  compared  with 
a Gaussian.  Whereas  in  class  G psfs  the  severity  of  the  blur  is  determined  by  how 
rapidly  the  optical  transfer  function  |/?.(fl)|  decreases  as  0 f 250.  in  defocus  blurs, 
severity  is  determined  by  the  number  of  zeroes1  in  |/?.(0)|  on  0 < 6 < 256.  Taking 
logarithms  of  absolute  values  in  Fig.  9B  emphasizes  the  difference  in  the  otf  signatures 
of  the  two  types  of  blur. 

Defocus  otfs  can  be  generated  in  Fourier  space  by  using  (18)  and  letting  the 
frequencies  £,r/  be  integers  with  —256  < £,77  < 250.  Specifying  R uniquely  determines 
the  defocus  otf.  Synthetically  defocused  images  can  be  generated  via  FFT  algorithms 


xThe  first,  five  positive  zeroes  of  Ji(x)  are  3.83171,  7.01559,  10.17347,  13.32369,  and  16.47063. 
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FOURIER  TRANSFORMS  OF  GAUSSIAN  AND  DEFOCUS  BLURS 


DISCRETE  FREQUENCY  ( INTEGER ) 
GAUSSIAN  ( DASHED ) DFFOCUS  {SOLID) 


B 


FOURIER  TRANSFORMS  OF  GAUSSIAN  AND  DE FOCUS  BLURS 


~-°  H 1 C 1 1 T 1 1 '~r 


-250  -1 50  -SO  SO  150  250 

DISCRETE  r ent  I EM)  ( INTEGER ) 

CAUSSIAN  (DASHED)  DEFOCUS  {SOLID) 


Fig.  9.  Fourier  domain  behavior  of  Gaussian  and  defocus  blurs.  A.  Comparison  of  e~a £ 
(dashed)  with  2Ji(R£)/(R£)  (solid),  on  —250  < £ < 250,  when  a = 0.001,  R = 0.08,  and  defocus  otf 
has  6 zeroes  in  0 < £ < 256.  B.  Comparison  of  corresponding  logarithms  of  absolute  values.  Zeroes 
in  defocus,  otf  appear  as  sharp  minima. 


Fig.  10.  Deblurring  with,  incorrect  defocus  parameter  and  honeycomb  artifacts.  A.  Noiseless 
defocused  Sydney  image  when  R = 0.1.  Blurred  image  computed  and  stored  in  64-bit  precision.  B. 
Successful  deblurring  of  image  A using  correct,  defocus  value  R = 0.1,  and  modified  SECB  procedure 
with  q(^,ri)  as  in  (34),  s — 0.001  and  1\  = 500.  C.  Deblurring  with  incorrect,  value  R = 0.11  produces 
‘ honeycomb ’ artifacts  when  A = 0.5,  despite  absence  of  noise  in  blurred  image.  With,  larger  errors 
in  R and/or  higher  values  of  I\  , artifacts  become  more  severe , eventually  obscuring  image  entirely. 


and  multiplication  with  2J\  (RQ)/(R0).  Deblurring  such  defocused  images  when  the 
defocus  parameter  R is  known,  can  be  done  using  the  modified  SECB  procedure 
discussed  in  the  Appendix.  This  is  a fast  FFT-based  direct  method. 

It  should  be  noted  that  the  defocus  variable  R used  in  this  paper  is  a Fourier 
domain  variable  that  enters  all  calculations  only  as  the  argument  of  the  Bessel  function 
J\  in  (18).  Our  numerical  values  for  R are  not  comparable  to  those  used  in  some 
studies  where  R denotes  the  radius  or  diameter  of  the  defocus  circle,  and  is  measured 
in  pixels.  In  comparing  the  severity  of  the  blurs  discussed  here  with  those  in  other 
studies,  the  relevant  measure  should  be  the  number  of  zeroes  of  h(9 ) on  0 < 0 < 25G. 

Our  first  example  is  the  defocused  Sydney  image  with  R — 0.1  in  Fig.  10A.  This 
was  calculated  and  stored  in  64-bit  precision  and  may  be  assumed  noiseless.  With 
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Fig.  11.  Influence  of  noise  on  defocus  parameter  identification  from,  blurred  image  trace.  A. 
Trace  of  defocused  Mariner  10  image  with  R = 0.08,  computed  and  stored  in  32-bit  precision  (solid 
line),  together  with  trace  of  defocus  otf  (dashed  line).  Sharp  local  minima  in  blurred  image  trace 
correctly  determine  defocus  parameter  R.  B.  Trace  of  defocused  Mariner  10  image  with  R — 0.08, 
stored  in  8-bit,  precision.  Pattern  of  sharp  local  minima,  now  seriously  affected  by  8-bit  rounding 
noise,  does  not  reliably  determine  R.  First  three  local  minima,  at  ( — 47,  £ = 90,  and  £ = 107, 
respectively  yield  R = 0.0815,  R = 0.0780,  and.  R = 0.095. 


the  otf  known  exactly,  excellent  deblurring  is  accomplished  in  Fig.  10B,  using  the 
modified  SEC'B  procedure  in  (33)  with  s = 0.001,  K = 500,  and  <y(£,  ?;)  as  in  (34). 
In  particular,  the  umbrellas  and  people  on  the  Opera  House  deck,  and  the  logos  on 
the  sailboats,  are  nicely  resolved  at  this  high  value  of  I\.  When  the  blurred  image 
contains  noise,  lower  values  of  I\  must  be  used,  resulting  in  lower  resolution. 

Of  great  significance  to  the  present  paper  is  the  image  in  Fig.  10C  that  results 
when  the  same  procedure  is  used  to  deblur  Fig.  10A,  but  with  the  incorrect  value 
R = 0.11.  The  honeycomb  artifacts  obscuring  the  image  in  Fig.  10C  are  not  typical 
noise  artifacts,  as  the  blurred  image  is  noiseless  and  I\  was  reduced  to  the  value 
A = 0.5.  Rather,  this  phenomenon  may  be  understood  as  follows.  The  SECB  formula 
in  (33)  may  be  written  as 


(19)  fHZ,v)  = k(tiri)9(Z,V), 


MCh) 

\kt,v)\2  + K~2\l  -qs{Z,T))\2' 


where  /c(£,  77)  is  the  Fourier  space  representation  of  the  regularized  inverse  to  the 
defocus  blur  operator.  The  blurred  image  |r)(C'//)|  Is  zero  or  *ias  very  small  values  at 
the  points  in  the  (£,  ry)  plane  where  | /?.(£,  ry)|  is  zero  or  has  very  small  values.  This  set  of 
points,  A0'1 , is  determined  by  the  value  R = 0.1.  With  that  value  of  R,  the  regularized 
inverse  |fc(£,ry)|  develops  appropriate  spikes  on  the  point  set  A01,  so  that  the  correct 
deblurred  image  results  upon  multiplication  of  A:(£,h)  with  <?(£,? y).  However,  with 
R = 0.11,  |A:(C  '/)i  develops  its  spikes  on  the  dislocated  point  set  A011,  where  the 
blurred  image  \g(f,r])\  need  not  always  be  small  or  zero.  Inevitably,  there  is  false 
overamplification  of  power  at  selected  frequencies  in  </(£,//)•  This  ‘resonance’  effect 
produces  the  artifacts  in  Fig.  10C.  With  larger  errors  in  R , and/or  larger  values  of  K, 
these  perturbations  change  character  and  become  more  intense,  eventually  obscuring 
the  image  entirely. 
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Several  techniques  for  blind  deconvolution  of  defocus  and  motion  blurs  seek  to  ex- 
ploit the  regular  patterns  of  zero  crossings  that  occur  in  the  optical  transfer  functions 
describing  such  blurs.  The  main  idea  is  illustrated  in  Fig.  11.  Using  B = 0.08,  the 
sharp  Mariner  10  image  in  Fig.  IB  was  synthetically  blurred  by  Fourier  multiplication 
with  (18),  and  the  resulting  image  was  stored  in  32-bit  precision.  At  this  low  level  of 
noise,  the  trace  of  log  |</*(£,0)|  in  Fig.  11A  (solid  curve)  displays  a regular  pattern 
of  sharp  minima,  at  the  same  locations  as  those  found  in  log  \2.J\  (0.08  #)/(0.08  D)\ 
(dashed  curve).  The  defocus  parameter  B can  easily  be  determined  from  blurred 
image  data  alone  in  that  case.  However,  when  the  blurred  image  is  stored  in  8-bit 
precision,  as  in  Fig.  11B,  log|y*(£,0)|  is  seriously  contaminated  by  8-bit,  round- 

ing noise,  and  the  regular  pattern  of  sharp  minima  is  no  longer  evident.  The  first 
three  local  minima  in  Fig.  11B.  at  £ — 47,  £ = 90,  and  £ = 107,  respectively  yield 
B = 0.0815,  B = 0.0780,  and  B — 0.095.  At  higher  noise  levels,  reliable  determination 
of  B from  blurred  image  data  alone  is  generally  not  possible. 

Using  additional  a-priori  information,  a more  sophisticated  and  computationally 
intensive  iterative  procedure  is  developed  in  [13],  [12],  The  method  is  based  on  view- 
ing the  blurred  image  as  a noisy  observation  of  an  autoregressive  moving  average 
(ARMA)  Markov  random  held.  Blur  detection  and  image  reconstruction  are  reformu- 
lated as  an  ARMA  parameter  identification  problem.  A maximum  likelihood  principle 
characterizes  the  desired  parameter  values  as  those  values  having  most  likely  resulted 
in  the  observed  image.  The  expectation-maximization  (E-M)  algorithm  is  then  used 
to  obtain  the  maximum  likelihood  solution.  The  procedure  requires  prior  estimates 
of  image  and  point  spread  function  support  sizes,  as  well  as  preprocessing  of  the 
boundaries  in  the  blurred  image.  The  authors  stress  that  the  E-M  algorithm  is  only 
guaranteed  to  converge  to  a stationary  point,  and  often  converges  to  one  of  several 
local  minima  rather  than  to  the  desired  global  minimum.  This  results  in  erroneous 
parameter  estimates  and  a failed  reconstruction.  The  convergence  point  is  strongly 
dependent  on  the  initial  guesses  for  parameter  values.  Repeated  trials  with  different 
initial  values  are  usually  required.  Good  a-priori  information  about  the  ideal  image 
and  the  blurring  kernel  is  essential  for  useful  results. 

The  method  described  below  is  a fast  FFT-based  direct  method.  At  the  levels  of 
noise  and  blur  intensities  where  it  produces  useful  reconstructions,  the  method  may  be 
used  in  two  different  ways.  It  may  be  used  as  a stand-alone  direct  blind  deconvolution 
technique,  or  it  may  be  used  in  conjunction  with  the  above  maximum  likelihood 
procedure,  or  some  other  iterative  approach,  to  provide  good  initial  values  and  other 
pertinent  information  necessary  for  convergence  to  the  desired  global  minimum. 

5.  A modified  BEAK  method  for  defocus  blur.  This  procedure  assumes 
knowledge  that  the  blur  is  a defocus  blur  and  requires  a substitute  image.  The  idea 
is  similar  to  that  in  section  3,  except  that  logarithms  are  not  used.  We  again  begin 
with  the  basic  relation 

(20)  g{x,y)  = h(x,y)  <g>  fe{x,y ) + n(x,y), 
where  the  noise  n(x,y)  is  assumed  to  satisfy 

(21)  I \n(x,y)\dxdy  <§C  f fe(x,y)dxdy  = cr  > 0, 

JR2  .1 R2 

so  that 


n*(f,»7)|  < 1- 
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Fig.  12.  Blind  deblurring  of  spacecraft  image  using  BEAK  psf  detection  method.  A.  BEAK 
detection  of  approximate  optical  transfer  function  in  8-bit  defocused  Mariner  10  image  (left),  us- 
ing gross  behavior  in  sharp  Mariner  5 image  (right,).  B.  Least  squares  fit,  of  trace  in  A with, 
2\.J\(Rf)/(Rf)\  on  —250  < f < 250,  returns  R = 0.081436,  not,  far  from,  exact  value  R = 0.08.  C.  De- 
focused Mariner  10  image.  D.  Successful  deblurring  of  image  C using  detected  value  R = 0.081436, 
and  modified  SECB  procedure  with.  q(fi,ri)  as  in  (3f),  s = 0.001  and  I\  = 0.5.  Faint  honeycomb 
artifacts  visible  against  dark  background. 


With  /?  > 0 and  ft  = yj £2  + rf , let  = {(£,  ?/)  | 6 < cc}  be  a neighborhood  of  the 
origin  where 

(23)  2 \J1(R9)/(R6)\  \.fe*(tv)\  » |n*(^»/)|. 

Such  an  Q exists  from  (22),  since  the  left  hand  side  of  (23)  is  unity  at  the  origin.  For 
(£,  //)  € we  have 

(24)  \g*(^v)\  « \h(Z,T])\  \fe*(tv)\  = 2 |Ji(i&9)/(ia9)|  \fe*(H,v)\. 

Hence,  for  |£|  < uj, 

(25)  \9*(Z,0)\  «2  |J1(i?0/(^)l  |//(£,0)|. 

We  now  replace  the  unknown  \ fe  (£,0)1  in  (25)  with  el’s^'1 , where  vs(f)  = — o|£|f'  is  the 
least  squares  approximation  to  log |/.s  (£,0)|,  and  fs(x,y)  is  the  substitute  image  for 
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Fig.  13.  Blind  deblurring  of  noisy  face  image  using  BEAK  psf  detection  method.  A.  BEAK 
detection  of  approximate  optical  transfer  function  in  noisy  defocused  Marilyn  Monroe  image  (left), 
using  gross  behavior  in  sharp  Ingrid  Bergman  image  (right.).  Exact  otf  has  9 zeroes  on  0 < 4 < 
256,  and  1%  noise  was  added  to  the  8-bit  blurred  image.  B.  Least  squares  fit,  of  trace  in  A with 
2\Ji(Rf)/(Rf)\  on  —150  < £ < 150.  returns  R — 0.118816,  not  far  from  exact  value  R = 0.12.  C. 
Noisy  defocused  Marilyn  Monroe  image.  D.  Successful  deblurring  of  image  C using  detected  value 
R — 0.118816,  and  modified  SECB  procedure  with  q(f,ii)  as  in  (34),  s = 0.001  and  K = 0.25. 
Visible  noise  artifacts. 


fe(x,  y )■  The  detection  procedure  is  then  the  following.  With  positive  a and  b given 
a-priori,  find  a positive  number  R such  that  2 |Ji(i?£)/(i?£)|  best  fits  |t/*(£,0)|  eal^ 
on  |£|  < co.  Note  that  the  expression 

(26)  u{0  = |$*(£,0)|  ea|fl\  |£l<^ 

which  involves  gross  behavior  in  a substitute  image,  is  a very  rough  approximation 
to  |fi(£,  0)|,  the  absolute  value  of  the  actual  defocus  optical  transfer  function  on  the 
line  i]  = 0.  Note  also  that  use  of  logarithms  would  require  the  fitting  of  log{«(£)} 
with  log{2  |Ji(-R£)/(-R£)|}-  The  sharp  spikes  in  the  latter  function  are  not  helpful  in 
the  fitting  algorithm.  Finally,  fitting  u(£)  produces  quite  different  plots  than  was  the 
case  in  Fig.  6A.  In  particular,  there  is  no  beak  at  £ = 0 in  this  version  of  the  BEAK 
method. 

We  illustrate  this  procedure  with  two  examples.  In  Figure  12  the  Mariner  10 
spacecraft  image  in  Fig.  IB  was  blurred  by  Fourier  multiplication  with  (18)  using 


16 


ALFRED  S.  CARASSO 


R = 0.08.  As  noted  in  Fig.  11  A,  the  otf  has  C zeroes  on  0 < £ < 256.  The  resulting 
binned  image,  rounded  to  8-bit.s,  has  the  trace  shown  in  Fig.  11B.  With  the  Mariner 
5 image  in  Fig.  1A  as  the  substitute  image,  we  have  a = 2.81,  b = 0.190.  We  now 
form  u(£)  in  (26).  The  plot  of  this  approximate  otf  is  shown  in  Fig.  12A.  Notice 
that,  the  maximum  value  in  u(£)  is  almost  5,  whereas  the  correct  value  in  the  true 
otf  is  1.  In  Fig.  12B  we  best  fit  u(£)  with  2 |Ji  (R£)/(R£)\  on  |£|  < 250.  This 
returns  7?  = 0.081436,  not  far  from  the  exact  value  R = 0.08.  Little  change  in  the 
returned  value  of  R results  when  the  fitting  procedure  is  used  on  the  smaller  interval 
|£|  < 100.  The  vertical  scale  in  Fig.  12B  is  restricted  to  a.  maximum  of  2 so  as  to 
better  display  the  resulting  fit  (solid  curve).  Using  the  modified  SECB  procedure 
with  <y(^,  //)  as  in  (34),  s — 0.001,  I\  = 0.5,  and  this  detected  value  of  /?.  produces  a. 
quite  good  reconstruction.  Faint  honeycomb  artifacts  are  in  fact  visible  against  the 
dark  background. 

Our  second  example  is  the  Marilyn  Monroe  image  in  Fig.  3A,  defocused  with  R = 
0.12.  In  that  case  the  otf  has  9 zeroes  on  0 < £ < 256.  The  blurred  image  was  rounded 
to  8-bits,  and  this  was  followed  by  the  addition  of  1%  uniformly  distributed  random 
noise,  i.e.,  each  8-bit  pixel  value  g{x,y ) was  replaced  by  {1  + 0.01  n(x,y)}g(x,y), 
where  n(x,  y)  is  a random  number  drawn  from  a uniform  distribution  in  the  range 
— 1 , 1] . The  Ingrid  Bergman  image  in  Fig.  3B  is  now  used  as  the  substitute  image. 
With  a = 3.11,  b = 0.183,  we  form  it(£)  in  (26).  This  approximate  otf,  displayed 
in  Fig.  13A  on  |<^|  < 150,  has  a maximum  value  of  about  3.  A least  squares  fit,  to 
these  data  with  2 | J\(R£)/ (R£)\  on  |£|  < 150  is  shown  in  Fig.  13B.  This  returns 
R = 0.118816,  not  far  from  the  exact  value  R = 0.12.  With  this  detected  value  of  R, 
the  modified  SECB  procedure  with  g(£,  //)  as  in  (34),  s = 0.001,  and  I\  = 0.25,  results 
in  a high  quality  reconstruction.  Noise  artifacts  are  clearly  visible  in  Fig.  13D. 

Similar  behavior  is  found  in  many  other  examples  of  moderately  defocused  images 
at,  low  levels  of  noise,  provided  a good  substitute  image  is  used.  However,  there  are  also 
many  cases  where  the  BEAK  method  does  not  provide  sufficiently  accurate  values  for 
R,  and  further  analysis  is  required  to  obtain  a useful  restoration.  The  next  example 
points  up  the  limitations  of  our  procedure,  and  signals  some  essential  difficulties  that 
attend  blind  deconvolution  by  any  method. 

6.  Failure  in  the  BEAK  method.  The  Mariner  10  experiment  in  Figure  12 
was  repeated  using  R = 0.25  when  the  defocus  otf  lias  20  zeroes  on  the  interval 
0 < £ < 256.  As  shown  in  Fig.  14A,  this  is  a,  severe  blur.  As  before,  the  blurred 
image  was  rounded  to  8-bits  and  the  Mariner  5 image  was  used  as  the  substitute 
image.  A least  squares  fit,  of  the  resulting  approximate  otf  with  2 | Ji{R.^) / (R,^)\  on 
|£|  < 50  returns  R = 0.22236.  Apparently,  this  is  insufficiently  close  to  the  correct, 
value  R = 0.25.  Restoration  with  that  detected  value  of  R produces  severe  honeycomb 
artifacts  that  obscure  the  image,  as  shown  in  Fig.  14B. 

While  the  presence  of  8-bit,  noise  in  the  blurred  image  is  not,  helpful,  noise  is  not 
the  primary  cause  of  failure  in  this  case.  Indeed,  if  the  R = 0.25  defocused  image  is 
computed  and  stored  in  64-bit,  precision,  and  the  above  procedure  is  reapplied  to  this 
noiseless  blurred  image,  the  maximum  detected  value  of  /?,  over  a number  of  trial 
subintervals  |£|  < u < 256,  is  R — 0.22468.  Deblurring  the  noiseless  image  using  that 
detected  value  of  R produces  results  very  similar  to  Fig.  14B,  even  when  A is  reduced 
to  the  value  I\  =0.1. 

Some  insight  into  the  reason  for  this  failure  is  provided  in  Figure  15.  In  Fig. 
15A,  which  is  analogous  to  Fig.  11A,  the  trace  of  log  \ge*(£i  0)|  (solid  line)  in  the 
64-bit  defocused  Mariner  10  image,  is  shown  on  the  interval  |£|  < 150,  together  with 
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Fig.  14.  Failure  of  BEAK  method  in  severely  blurred  image.  A.  8-bit.  defocused  Marnier  10 
image  when  R — 0.25  and  off  has  20  zeroes  on  0 < £ < 256.  BEAK  detection  of  approximate 
transfer  function  using  sharp  Mariner  5 image,  and  subsequent,  least  squares  fit  on  —50  < £ < 50, 
returns  R = 0.22236,  insufficiently  close  to  correct  value  R ~ 0.25.  B.  Deblurring  of  image  A using 
detected  value  R = 0.22236,  and  modified  SECB  procedure  with  q(f,ri)  as  in  (34),  s = 0.001  and 
K = 0.25,  produces  severe  honeycomb  artifacts  that  obscure  image. 

A B 


64-BIT  DEFOCUSED  MARINER  10  IMAGE  WHEN  R-t)  25 


64-BIT  DEPOCDSED  MARINER  10  IMAGE  WREN  R-0  25 


DISCRETE  FREQUENCY  {INTEGER) 

DEFOCIIS  PSF  {DASHED  I BLURRED  IMAGE  (SOLID) 


DISCRETE  FREQUENCY  (INTEGER) 

DEPOCUS  PSF  (DASHED)  BLURRED  IMACE  (SOLID) 


Fig.  15.  Misleading  patterns  in  noiseless  blurred  image  data  give  false  reading  for  R.  A. 
log|ffe*(£,0)|  in  64-bit  defocused  Mariner  10  image  when  R = 0.25  (solid),  and  log  |/i(£,  0)|  (dashed). 
Despite  absence  of  noise,  sharp  local  minima  in  blurred  image  trace  are  not  well-correlated  with 
zeroes  of  .J i(0.25£).  B.  Identifying  first  3 sharpest  mamma  in  noiseless  blurred  data  with  first  3 
positive  zeroes  of  Ji(x)  implies  R ~ 0.15,  grossly  underestimating  true  value  R = 0.25. 


the  logarithm  of  the  defocus  otf  when  Ft  = 0.25  (dashed  line).  Evidently,  there  are 
many  more  local  minima  in  the  defocus  otf  in  Fig.  15 A,  as  compared  to  Fig.  11  A, 
even  though  the  ^-interval  in  Fig.  15A  has  been  reduced  to  improve  visibility.  This 
produces  a crowding  effect.  In  addition,  the  defocus  otf  and  blurred  image  data  are 
only  available  at  discrete  values  of  f.  For  a given  It.  the  zeroes  of  J\  (Rf  ) cannot 
always  be  well-captured  on  a preassigned  discrete  mesh.  A zero  of  -J\  that  does  not 
fall  on  a mesh  point  translates  into  a small  value  of  the  defocus  otf  at  the  nearest 
mesh  point,  if  it  is  not  missed  altogether.  And,  if  the  otf  value  at  that  mesh  point 
is  not  sufficiently  small,  the  resulting  local  minimum  in  the  logarithmic  defocus  trace 
will  not  be  particularly  sharp. 
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A 


B 


COMPARISON  OF  BEAL  DETECTED  OTF  KITH  EXACT  OTF 


COMPARISON  OF  FITTED  OTP  KITH  EXACT  OTF 


Fig.  16.  Analysis  of  failure  of  BEAK  method  in  severely  blurred  mage.  A.  Comparison  of 
approximate  otf  u(£)  (dashed  curve)  with  exact  otf  ueR)  (solid  curve)  on  |£|  < 30.  Minima  in 
u(f)  at  f = 15  and  £ = 28,  coincide  with  those  in  ue((,)  and  are  compatible  with  R ss  0.25.  B. 
Comparison,  of  least  squares  fit  to  u(f)  (dashed  curve)  with  ite(£)  (solid  curve)  on  |£|  < 50.  Minima 
in  least  squares  fit  dislocated  relative  to  true  values.  Fit  returns  R = 0.22236,  locating  first  minimum 
at  f ~ 17  rather  than  at  £ = 15. 

This  behavior  of  the  defocus  otf  at  preassigned  discrete  values  of  £ is  translated 
into  similar  behavior  of  the  blurred  image  trace  at  these  same  £ values.  The  blurred 
image  is  the  primary  information  in  blind  deconvolution,  and  the  discrete  mesh  was 
imposed  when  the  digitized  image  was  acquired.  Weak  local  minima  in  log|/i(£,0)| 
become  weak  local  minima  in  log  | ge*(£,  0)|.  As  a result,  the  crowded  minima  resulting 
from  zeroes  of  J\  no  longer  clearly  stand  out  among  other  local  minima  in  the  highly 
oscillatory  blurred  image  trace.  Hence,  despite  the  absence  of  noise  in  the  64-bit  solid 
curve  in  Fig.  15A,  there  is  no  convincing  regular  pattern  of  sharp  minima  analogous 
to  that  in  Fig.  11  A. 

This  point  is  more  clearly  made  in  Fig.  15B,  which  displays  log|h(£,0)|  and 
log  |</e*(£,  0)1  on  the  interval  0 < £ < 70.  Without  the  knowledge  provided  bv  the 
dashed  curve,  the  relatively  sharp  minima  in  the  solid  curve  at  £ = 28,  £ = 47  and 
£ = 66,  and  respectively  labeled  1,  2 and  3,  might  easily  be  construed  as  originating 
from  the  first  3 positive  zeroes  of  Ji(i?£).  Such  a scenario  yields  3 distinct  but  not 
seriously  incompatible  values  for  R,  namely,  R = 0.1368,  R = 0.1493  and  R = 0.1541. 
In  fact,  very  similar  considerations  are  used  to  initiate  iterative  blind  algorithms, 
and  any  one  of  these  3 values  of  R might  be  a plausible  initial  guess  in  the  ARMA 
maximum  likelihood  approach  in  [13],  [12],  for  example.  Convergence  to  the  desired 
global  minimum  from  such  initial  values  would  be  doubtful. 

Clearly,  any  detection  procedure  that  examines  the  solid  curve  in  Fig.  15A  and 

concludes  that  R = 0.25,  must  do  so  on  the  basis  of  quite  good  a-priori  information. 

~ * 

Very  helpful  but  unavailable  a-priori  knowledge  would  be  log  \ f€  (£,  0)|,  the  trace  in 
the  exact  sharp  Mariner  10  image,  since  subtracting  the  latter  from  the  solid  curve 
in  Fig.  15  A produces  log  | // ( £ , 0)|,  the  dashed  curve  in  Fig.  15A.  Equivalently,  we 
may  exponentiate  these  traces  and  form 

(27)  M£)  = l<?e*(£,0)|/|/e*(£,0)|  = |/?.(£.  0)  | . 

The  BEAK  method  seeks  to  emulate  the  process  leading  to  (27)  using  available  infor- 
mation. With  the  noisy  data  |</*(£,0)|  and  the  gross  behavior  us(£)  = — o|£|6  in  the 


BLIND  DECONVOLUTION  II 


19 


substitute  Mariner  5 image,  we  form  the  approximate  otf  »(£)  in  (26).  The  trace  of 
u(£)  on  |£|  < 30  is  shown  in  Fig.  IGA  (dashed  curve),  and  compared  to  the  exact  otf 
ue (£)  (solid  curve).  It  is  noteworthy  that  «(£)  has  local  minima  at  £ = 15  and  £ = 28. 
These  minima  coincide  with  the  first  two  minima  in  ■«,.(£)  and  are  compatible  with 
R ~ 0.25.  (The  first  positive  zero  in  ,/i(0.25£)  occurs  at  £ = 15.32684  which  is  not 
a mesh  point).  However,  the  least  squares  fit  to  u(£)  with  2 \Ji(R^)/(R^)\  does  not 
coincide  with  ue(£)  !!  Instead,  the  more  objective  examination  of  the  u(£)  data  by 
the  least  squares  algorithm  results  in  the  dashed  curve  in  Fig.  I6B,  whose  minima 
are  dislocated  relative  to  those  in  the  exact  otf  (solid  curve).  The  returned  value  for 
R is  0.22236,  and  the  fit  locates  the  first  minimum  at  £ = 17,  rather  than  at  £ = 15. 

7.  Possible  recovery  from  failure.  While  the  behavior  in  Fig.  14B  is  more 
common  in  severely  defocused  images,  similar  behavior  due  to  inaccurately  detected 
R values  can  occur  with  more  moderate  blurs.  There  are  three  avenues  that  can 
be  explored  to  improve  reconstruction  in  such  cases.  One  should  first  compare  the 
first  two  or  three  minima  in  the  fitted  curve  with  the  corresponding  minima  in  u(£), 
provided  the  latter  minima  are  clearly  defined.  Comparing  logarithms  of  both  curves 
is  often  helpful.  If  there  is  dislocation,  new  values  for  R corresponding  to  the  minima 
in  u(£)  can  be  tried.  A second  approach  assumes  the  BEAK  returned  value  Rs  to  be 
within  striking  distance  of  the  exact  value  Rc.  For  512  x 512  images,  20  trial  SECB 
restorations,  each  with  a different  value  of  A.  can  be  obtained  in  about  a minute  of 
cpu  time  on  current  desktop  workstations.  Starting  from  Rs , an  interactive  visual 
search  for  Re  is  therefore  quite  feasible.  In  that  search,  care  must  be  taken  to  use 
low  values  for  K in  the  SECB  method,  until  a value  of  R is  located  where  artifacts 
disappear.  Higher  values  of  A can  then  be  used  to  improve  resolution.  Finally,  R n 
can  be  used  as  the  initial  guess  in  an  iterative  blind  deconvolution  algorithm.  In  the 
case  of  Fig.  14B,  the  BEAK  returned  value  Rb  — 0.222  is  a more  useful  initial  guess 
than  are  the  values  R ph  0.15  obtained  by  inspecting  the  solid  curve  in  Fig.  15B. 

8.  APPENDIX.  A modified  SECB  procedure  for  psfs  not  in  G.  The 

SECB  method  applied  to  the  convolution  integral  equation  H f — g in  (3)  presupposes 
the  psf  h(x , y)  to  have  a non-vanishing  Fourier  transform  h(£,  g),  so  that  for  any  fixed 
s > 0,  /is(£,  y)  can  be  uniquely  defined  [7,  p.  555].  In  addition,  /?s(£,  y)  is  required  to  be 
the  Fourier  transform  of  a probability  density.  This  is  the  case  for  class  G point  spread 
functions  in  (4),  and  more  generally,  for  2-D  infinitely  divisible  probability  densities, 
[7].  The  integral  operator  Hs  is  then  defined  to  be  the  operator  of  convolution 
with  the  inverse  transform  of  hs(£,  y).  Solving  (3)  is  equivalent  to  solving  a time- 
reversed  diffusion  equation,  and  can  be  implemented  as  a continuation  problem  using 
the  operator  Hs.  Such  a continuation  approach  is  an  important  element  of  the  APEX 
method  discussed  in  [6].  The  SECB  constrained  solution  /I  to  the  ill-posed  equation 
Hf  = g is  defined  by 


(28)  fHx,y)  = Arg  ( min  (||  Hf  - g 

l feL 2 


+A  - II  H*f  - f 


where  the  positive  constants  A , s are  previously  chosen  regularization  parameters,  as 
discussed  in  [2,  3,  4],  This  minimization  problem  has  a closed  form  solution  in  Fourier 
space.  We  have,  with  ~z  denoting  the  complex  conjugate  of  z. 


MCh)fl(Ch) 

\h{£,  y)\2  + A’-2|l  — hs(£,y)\2' 


(29) 
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Thr  functional  on  the  right,  of  (28)  embodies  the  a-priori  information  that  with  suitable 
fixed  K,  s > 0,  the  correct  solution  f°(x,y)  satisfies 

(30)  \\  H f°  — g \\<  e,  ||  Hsf°  - f°  ||<  AT, 

where  e > 0 is  an  L1  bound  for  the  noise  in  g(x,y).  Rigorous  error  bounds  for  SECB 
regularization  are  developed  in  [3.  4].  Applications  of  the  SECB  constraint  to  ill-posed 
problems  other  than  image  deblurring  are  discussed  in  [5]. 

The  above  diffusion  formalism  is  not  applicable  to  out  of  focus  blurs  since  the 
continuous  real-valued  function  //.(£,•//)  hr  (18)  has  infinitely  many  sign  changes.  The 
operator  Hs  is  not  well-defined,  and  solving  (18)  is  no  longer  equivalent  to  solving 
a diffusion  equation  backwards  in  time.  For  deconvolution  problems  Hf  = g where 
//(£,//)  is  not  an  infinitely  divisible  characteristic  function,  we  may  proceed  as  follows. 
Choose  i/(£, rj)  = exp{—a(f2  + r/~ ) ' } with  fixed  o > 0 and  fixed  0 < (3  < 1,  and  let 
Qs  denote  convolution  with  the  inverse  transform  of  qs{f,q),  0 < s < 1.  The  correct 
solution  f°(xpy)  of  Hf  = g again  satisfies 

(31)  ||  Hf°  -g\\<e,  ||  Qsf°  - f°  || < AT, 

with  suitably  chosen  A",  s > 0.  Accordingly,  we  obtain  a constrained  solution  /'  of 
H f = g by  means  of 

(32)  p(x,y)  = Ary  | min  (||  Hf  - g ||2  +K~2  ||  Qs f - / ||2)|  , 
leading  to  the  Fourier  domain  solution 


(33) 


fHZ,v) 


h(Ch)g(Ch) 

g)\2  + A 1 1 - qs(f,r])\2 


The  error  analysis  in  [3]  must  be  modified  in  order  to  be  applicable  in  (32).  For  a 
given  fixed  rj),  use  of  (33)  on  a blurred  image  g(x,  y)  requires  knowledge  of  the 
regularization  parameters  A',  s.  These  parameters  will  depend  on  the  choice  of  Q 
and  represent  a-priori  information  about  the  exact  solution,  some  form  of  which  is 
always  necessary  in  the  solution  of  ill-posed  problems.  In  practice,  if  the  unknown 
sharp  image  is  an  easily  recognizable  object,  we  may  fix  a value  of  s in  the  range 
0.U01  < s < 0.01,  and  adjust  K interactively  in  (33)  so  as  to  achieve  optimal  results. 
Higher  levels  of  noise  dictate  smaller  values  of  A",  and  vice-versa.  Begining  with  a 
small  value  of  A’,  increasing  I\  increases  sharpness  in  the  restored  image,  until  a 
threshold  value  is  reached.  Further  increase  in  I\  brings  out  noise  which  eventually 
obscures  the  image.  In  the  defocus  examples  in  Sections  4-0,  q(f,  r/)  was  chosen  to  be 


(34) 


q(^v) 


e-0.075(eW)1/2 


with  integer  £,  q lying  between  —256  and  250. 
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